%IEEE14 test
Edgetable = table();
Edgetable.EndNodes = [
                      1 2;
                      1 5;
                      2 3;
                      2 4;
                      2 5;
                      4 3;
                      4 7;
                      4 9;
                      5 4;
                      5 6;
                      6 11;
                      6 12;
                      6 13;
                      7 9;
                      8 7;
                      9 10;
                      9 14;
                      11 10;
                      12 13;
                      13 14;
                      ];

Edgetable.P2 = [
               77.95;
               42.05;
               32.98;
               35.85;
               27.41;
               1.22;
               12.42;
               11.33;
               36.92;
               24.95;
               7.43;
               7.71;
               17.61;
               32.42;
               20.0;
               5.07;
               9.18;
               3.93;
               1.61;
               5.72
               ] * 1000;

Edgetable.P1 = Edgetable.P2;

Pg = [
      1 120.0 875e-6;
      2 40.0 525e-6;
      3 60.0 0;
      6 19.0 520e-6;
      8 20.0 0;
      ];
Pg(:, 2) = Pg(:, 2) * 1000;

Pd = [
      2 21.7;
      3 94.2;
      4 47.8;
      5 7.6;
      6 11.2;
      9 29.5;
      10 9.0;
      11 3.5;
      12 6.1;
      13 13.5;
      14 14.9;
      ];
Pd(:, 2) = Pd(:, 2) * 1000;

% Pinj = zeros(14, 1);
% Pinj(Pg(:, 1)) = Pg(:, 2);
% Pinj(Pd(:, 1)) = Pinj(Pd(:, 1)) - Pd(:, 2);
% Nodetable = table(Pinj);

Nodetable = table(zeros(14, 1), zeros(14, 1), 'VariableNames', {'Pg', 'Pd'});
Nodetable.Pd(Pd(:, 1)) = Pd(:, 2);
Nodetable.Pg(Pg(:, 1)) = Pg(:, 2);

%潮流分布图
PowerGraph = digraph(Edgetable, Nodetable);

%分布矩阵，路径矩阵，节点有功通量
[H, D, Psum] = power_distribution(PowerGraph);
PowerGraph.Nodes.Psum = Psum;

%节点碳流通量
PowerGraph.Nodes.Csum = (D.' * (Pg(:, 3) .* Pg(:, 2)));

%支路碳流率
Ce2 = zeros(numedges(PowerGraph), 1);

for k = 1:numedges(PowerGraph)
    s = PowerGraph.Edges.EndNodes(k, [1 2]);
    Ce2(k) = PowerGraph.Nodes.Csum(s(1)) * H(s(1), s(2));
end

PowerGraph.Edges.Ce2 = Ce2;

%支路碳流密度
PowerGraph.Edges.Cr = PowerGraph.Edges.Ce2 ./ PowerGraph.Edges.P2;

%负荷碳流率
PowerGraph.Nodes.Cd = PowerGraph.Nodes.Pd ./ PowerGraph.Nodes.Psum .* PowerGraph.Nodes.Csum;

%机组至负荷碳流率，不在表里的说明不存在通路
s=size(Pg,1)*size(Pd,1);
Cg2Cd=struct();
Cg2Cd.G=zeros(s,1);
Cg2Cd.L=zeros(s,1);
Cg2Cd.P=zeros(s,1);
Cg2Cd.C=zeros(s,1);
nidx=0;
for kg=1:size(Pg,1)
    for kd=1:size(Pd,1)
        d=D(kg,Pd(kd,1));
        if d>0
            nidx=nidx+1;
            Cg2Cd.G(nidx)=Pg(kg,1);
            Cg2Cd.L(nidx)=Pd(kd,1);
            Cg2Cd.P(nidx)=Pg(kg,2)*d*Pd(kd,2)/PowerGraph.Nodes.Psum(Pd(kd,1));
            Cg2Cd.C(nidx)=Cg2Cd.P(nidx)*Pg(kg,3);
        end
    end
end
Cg2Cd.G=Cg2Cd.G(1:nidx);
Cg2Cd.L=Cg2Cd.L(1:nidx);
Cg2Cd.P=Cg2Cd.P(1:nidx);
Cg2Cd.C=Cg2Cd.C(1:nidx);
Cg2Cd = struct2table(Cg2Cd);

%% 画图

%分布因子矩阵
f=figure('Units','centimeters','Position',[6,10,10,7]);
h1=heatmap(H);
Eula=[052 073 100;177 217 229;172 178 188;246 246 246;071 094 165]./255;
cmap2=interp1((1:3)',Eula([4 2 1],:),linspace(1,3,400));
colormap(h1,cmap2);
g=gca;
set(g,'FontSize',10,'fontname','Times new Roman'); 
set(g, 'Units', 'normalized', 'OuterPosition',[-0.06,-0.04,1.10,1.08]);


%机组-节点矩阵
f=figure('Units','centimeters','Position',[6,10,10,7]);
h2=heatmap(1:14,compose('G%d',1:5),D);
colormap(h2,cmap2);
g=gca;
set(g,'FontSize',10,'fontname','Times new Roman'); 
set(g, 'Units', 'normalized', 'OuterPosition',[-0.06,-0.04,1.10,1.08]);

%桑基图
figure('Name','power distribution','Units','normalized','Position',[.05,.2,.5,.56])
SK=SSankey(compose('G%d',Cg2Cd.G),compose('%d',Cg2Cd.L),table2cell(Cg2Cd(:,3)));
SK.draw();
figure('Name','carbon distribution','Units','normalized','Position',[.05,.2,.5,.56])
idx=Cg2Cd{:,4}>0;
SK=SSankey(compose('G%d',Cg2Cd.G(idx)),compose('%d',Cg2Cd.L(idx)),table2cell(Cg2Cd(idx,4)));
SK.draw();
